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Abstract 

In this paper we present a multigrid approach to accelerate the convergence of the iterative method 
proposed in [T3] to solve the Poisson equation in arbitrary domain fl, identified by a level set function ip, 
Q — {x £M.'^: ip{x) < O}, and mixed boundary conditions. The discretization is based on finite difference 
scheme and ghost-ccU method. This multigrid strategy can be applied also to more general problems 
where a non-eliminated boundary condition approach is used. Arbitrary domain make the definition of 
the restriction operator for boundary conditions hard to find. A suitable restriction operator is provided 
in this work, together with a proper treatment of the boundary smoothing, in order to avoid degradation 
of the convergence factor of the multigrid due to boundary effects. Several numerical tests confirm the 
good convergence property of the new method. 



Introduction 

Multigrid technique is one of the most efficient strategy to solve a class of partial differential equations, using 
a hierarchy of discretizations. It accelerates the convergence of an existing iterative method, which otherwise 
slowly converges toward the solution of the discrete problem, due to the bad convergence rate for the low 
frequency components of the error. The idea of multigrid method is to solve such low frequency component 
in a coarser grid. An introduction to multigrid can be found, for example, in [4j, while more advanced 
textbook on the subject are, for example, |41L I21j. Most iterative schemes to solve Elliptic equations can 
be speeded up by a multigrid technique. 

Elliptic equation in arbitrary domain (possibly with moving boundary) is central to many applications, 
such as diffusion phenomena, fluid dynamics, charge transport in semiconductors, crystal growth, elec- 
tromagnetism and many others. The wide range of applications may require different kind of boundary 
conditions. Let us look for instance at the temperature distribution in a medium of arbitrary shape satisfy- 
ing stationary heat equation: we may have Dirichlet (the temperature is fixed at the boundary), Neumann 
(heat flux is prescribed), or mixed boundary conditions (namely different boundary conditions on different 
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part of the boundary). More general Robin boundary conditions may also be prescribed, as in Stefan-type 
problem, in which a combination of temperature and heat flux is prescribed at the boundary (e.g. see [201 [S]). 
An application we have in mind is to fluid dynamics: the aim is to model the motion of an incompress- 
ible fluid contained in a tank of arbitrary shape. The problem is modeled by incompressible Navier-Stokes 
equations, which are solved by projection method of Chorin [101 El]. This leads to an elliptic equation for 
the pressure, obtained enforcing the incompressibility condition. This pressure equation requires Dirichlet 
condition on the free surface of the fluid and Neumann condition on the rigid walls. The pressure equation 
is the bottleneck of the whole method and therefore requires an efficient solver. 

Several techniques have been developed to solve Elliptic equation on an arbitrary domain. Finite Element 
Methods use a mesh triangulation to capture the boundary, such as in [34^ l35l [30l |37] . However, in presence 
of moving boundary, a grid re-meshing is needed at each time step, which makes the method expensive. 
Furthermore, for a complex geometry, generation of a good mesh is a non trivial task that may require a 
considerable amount of work [28j . For this reason a Cartesian grid method is preferred together with a 
level-set approach to keep track of the boundary at each time step. Level-set methods have been introduced 
to implicitly define a domain and its boundary, in order to simple handle complex topological changes of 
moving boundary such as merging and breaking up. Several papers and books exist in the literature about 
level-set method: [iQl [311 [36l [T71 [25] are just some examples. 

Since the boundary may be not aligned with the grid, a special treatment is needed. The simplest 
method makes use of the Shortley-Weller discretization [38], that discretizes the Laplacian operator with 
usual central difference away from the boundary and makes use of a non symmetric stencil in the interior 
points of the domain close to the boundary. While this discretization provides a simple second order 
method for Dirichlet conditions, it cannot be immediately applied in presence of Neumann conditions. In 
fact, Shortley-Weller discretization |38j for Neumann conditions requires that the value of the numerical 
solution is suitably reconstructed at the intersection between the grid and the boundary by applying the 
boundary condition. This approach is adopted, for example, by Hackbusch in j22j to first order accuracy, 
and by other authors (see [2] and the references therein) to second order accuracy. However, the method 
proposed by Bramble in [2] for second order accuracy is quite involved and not recommendable for practical 
purposes. 

Another class of methods is based on cut-cell methods, obtained by a Finite Volume discretization which 
embeds the domain in a regular Cartesian grid |24] . Cells that are cut by the boundary requires a special 
treatment, such as cell-merging and rotated-cell, in order to avoid a too strict restriction of the time step 
dictated by the CFL condition (e.g. see [23l[9l [T2] ) . 

Other methods for Dirichlet condition are the Immersed Boundary Method, first proposed by Peskin 
in [33], and later developed by several other authors [271 112]) with a proper multigrid approach [Ij, and 
penalization methods [7]. 

In our method we will use a rather simple finite-difference ghost-cell technique, that adds extra grid 
points (ghost points) outside the domain in order to keep unchanged the symmetry of the stencil even for 
inside points close to the boundary. A detailed description of the method can be found in [H] . 

In ghost points the boundary conditions are enforced in order to close the discrete system. The ghost- 
cell method was first developed by Fedkiw in [18] , where a two-phase contact discontinuity was discretized 
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(Ghost Fluid Method). A second-order accurate method for Dirichlet conditions on regular Cartesian grid 
is proposed by Gibou et al. in [19] . The value at the ghost nodes is assigned by linear extrapolation, and 
the whole discretization leads to a symmetric linear system, easily solved by a preconditioned conjugate 
gradient method. A fourth order accurate method is also proposed in [20j. Other methods use a non-regular 
Cartesian grid, such as in |8], where Gibou et al. present finite difference schemes for solving the variable 
coefficient Poisson equation and heat equation on irregular domains with Dirichlet boundary conditions, 
using adaptive Cartesian grids. One efficient discretization based on cut-cell method to solve more general 
Robin conditions is proposed by Gibou et al. in [32j, which provides second order accuracy for the Poisson 
and heat equation and first order accuracy for Stefan- type problems. 

Most of the techniques listed above cannot be straightforwardly applied in the special case of mixed 
boundary conditions. For cut-cell based methods I32j . different boundary conditions cannot to be 
imposed on the same boundary edge of a cut cell. Simple efficient methods based on symmetric image of 
ghost points to solve mixed boundary condition problems provided with a multigrid algorithm have been 
recently developed in [6] and by Ma et al. in [29]. 

In our method [13], boundary conditions are neither eliminated from the discrete system (they are 
strongly coupled and their elimination is too hard to perform in more than one dimension) nor directly 
enforced (which leads to a non-convergent iterative method): they are relaxed together with the interior 
equations. This leads us to an iterative scheme for the set of all unknowns (internal points and ghost 
points), which is proved to converge, at least for first order accurate discretization. 

In this paper we provide a general multigrid technique to solve the discrete system coming from a 
continuous elliptic problem in case of non-eliminated boundary conditions. The smoothing procedure of the 
multigrid approach in the interior is Gauss-Seidel-like, while the iterations on the boundary are performed 
in order to provide smooth errors. 

Multigrid techniques for non-eliminated boundary conditions are well-studied in literature in the case 
of rectangular domain (as we can see in |21|, I41j). where a restriction operator is defined separately for 
the interior of the domain and for the boundary, and the restriction of the boundary is performed using a 
restriction operator of codimension 1, since ghost points are aligned with the Cartesian axis. In the case 
of arbitrary domain, ghost points have an irregular structure and we provide a reasonable definition of the 
restriction operator for the boundary conditions. The method proposed in this paper can be extended to 
the case of discontinuous coefficient: a preliminary result in one dimension can be found in [14J, while the 
two-dimensional case is in preparation. 

In this paper we also show that a proper treatment of the boundary iterations can improve the rate 
of convergence of the multigrid, making it closer to the convergence rate predicted by the Local Fourier 
Analysis for inside equations, as suggested by Brandt in [3]. The cost of this extra computational work is 
negligible, i.e. tends to zero as the dimension of the problem increases. A comparison with other kinds of 
treatment of the boundary condition smoothing procedure is carried out. 

The paper is divided in three sections. We start with the multigrid approach in the one dimensional 
case, described in Section 1, with a special treatment of the transfer operators. Most of this method can be 
extended to high dimension, treated in Section 2, but a special care has to be taken for transferring the defect 
of boundary conditions. In this section a level-set approach is also introduced. Section 3 provides a strategy 
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to improve the convergence factor making it closer to the one predicted by the Local Fourier Analysis for 
interior relaxations. Numerical evidence of such improvement is provided, together with a comparison with 
other boundary condition smoothers (Kaczmarz and Block relaxation). 
In all the paper, we mainly use the notation of [41j . 



1 One-dimensional case 

In this section we will follow the description of the method proposed in fl3|, which is second order accurate, 
and provides a multigrid approach to speed up the convergence of the iterative scheme. For one-dimensional 
case, the multigrid approach in arbitrary domain is a natural extension of the basic multigrid strategy that 
can be found in any good basic text about multigrid, such as |4H [H [21] . Although if we can eliminate the 
boundary conditions from the linear system obtained by discretizing the problem, we always want to treat 
the case of non-eliminated boundary conditions in order to straightforwardly extend the method to more 
than one dimension, where the elimination of the boundary conditions from the system is hard to perform. 



1.1 Model problem and relaxation scheme 

Let D = [—1, 1] be the computational domain, a and h constants such that — 1 < a < 6 < 1, and = [a, h\. 
Let > 1 be a fixed integer and h = 2/N the spatial step, let Dh = {—1 = xq < xi < . . . < xn = 1} be 
the set of equally spaced grid points, and 17^ = n the set of inside grid points. Consider the model 
problem 

-u" = f inn 

u{a) = Qa (1) 
u'{h) = gi,. 

Let / and r be such that xi < a < x^+i, Xr~i < b < Xr (see Figure [l]). We use a ghost-cell method 
to discretize the problem. In order to obtain an iterative method, we solve the associate time- dependent 
problem 



du{t, a) 
dt 



liD{9a- u{t,a)) (3) 



du{t,b) ( du{t,by. 

- I^Nlgb w— (4) 



dt \ dx 

u{0,x) = uq{x) in fl (5) 

and we look for a second order accurate steady state solution, which is the solution of the original problem. 

Let us begin to discretize ([2]) in ri/^ = . . . ,Xr-i}- We use central difference in space and forward 

Euler in time for ^ obtaining: 



^ ^(H ^ At _ ^ ^ ^ ^^^^^ i = l + l,...,r-l. (6) 
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Taking the maximum time step consented by CFL condition [15], i.e. At = h? /2, we obtain: 

ut^'^ = \{utl+u\':l+ef), ^ = l + l,...,r-l. (7) 

Note that if we discretize directly the first equation ([T]) using central difference for the Laplacian operator, 
and use Jacobi iterative scheme for such discretization, we obtain exactly ([T]). 



-1 a^^^'^ ^^^^ ^^'^ ^'''^b ^ 

Fig. 1: Discretization of the domain in ID 

To obtain second order accuracy, we have to discretize the spatial terms in (|3]), Q to second order, while 
first order time discretization can be used, because we are just interested at the accuracy as t — )• +oo. 

We then can use linear interpolation for u{t,a) in Q. Since in some application it is required second 
order accuracy of the gradient of the solution, we use quadratic interpolation instead linear interpolation, 
obtaining: 

nf-"-'^ = - f^oAt ((1 + ^Ol^;™^ + (1 + W - ^i)ut:l - (1 - - 5a) , (8) 

where i9/ = — a)/h, and quadratic interpolation of u in nodes for the du{t,b)/dx in Q, 

obtaining: 

_ „{m) _ f^N At f (rn) _ (m) / (m) _ (m) („) 



ur^'> = ur> - ^ [ui^i - + - 2^ri + j + j j+t^N^t g,. (9) 

The constants fio and are chosen in order to satisfy the CFL conditions, i.e. fioAt < 1 and njyAt/h < 
2/3 (see [H])- In numerical tests of Section [3| we choose /^uAt = 0.9 and fi^At = 0.9 • 2/i/(3). Since 
At = h'^/2, then = 1.8/h'^ and jjn = 3.6/(3/i). 

In summary, our second order accurate iterative method is described by Eqs. Q, Q and ([9]), with the 
choice of constants 

At = h^/2, = 1.8//^^ /iTV = 1.2//1. (10) 

1.2 Multigrid approach 

We call F/j the set of ghost points, i.e. = {xi,Xr}. Let Ih be a general subset of Dh. We introduce the 
linear space of grid functions over and we denote it S{Ih) = {"Wh: //^ — M}. For any Wh G S{Ih), we 



5 



pose = w/i(xi). Let ih € S{^lh) such that fj^ = f{xi). The iterative scheme ([7|, (|8|, ^ converges to the 
exact solution of the discretized system 



9a 



S'iv(uft) = 9b, 



(11) 
(12) 

(13) 



where A^: S{Qh U T^) — >■ S{Qh) is defined by: 



AhUh{xi 



''i-l 



2«' + <i 



/l2 



while g]j,g% ■ S{nh U Th) — )• M are the discrete versions of the boundary conditions: 
g^UH) = (1 + + (1 + ^0(1 - - (1 - ^0|<2, 



^-2 



''r-2 



h 



2u: 



r-l 



System ( 11 )-( 13 ) can be interpreted in general as a discrete system of a Poisson equation with non-eliminated 
boundary conditions. 

Let us consider an arbitrary grid function G S{Qh U Th) and let 



9a 

9b 



9a-9Uu') 
9b-9%in'') 



be the defects of (11), (12), (13) respectively. Because of the linearity of Ah, g^, g%, if we solve exactly the 
so-called residual problem 



9D{^h) = 9a 
9N{^h) = gb 



(14) 

(15) 
(16) 



in the unknown € S{U,h U F/j), then \ih = + e/j is the exact solution of the system (11), (12), (13). In 
the basic idea of multigrid one needs to solve the residual problem in a grid coarser than the original one. 
We can summarize the iterative scheme ([T]), ([s]), ^ as follows: 



u 



(m+l) 



[^t^\h,9a,gb^ 
■■ s{^h u Th) X s{vth) X m2 ^ siVLh u r^. 



(17) 
(18) 
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Note that the iterative scheme ([T]), ([s]), ^ is of a Jacobi kind. In order to provide a multigrid strategy, 
we just require that the iteration operator (18) has the smoothing property, i.e. after few iteration steps 



(17), the error becomes smooth (not necessarily small). Roughly speaking, the high-frequency components 



of the error reduce quickly. We call smoothers any operator (18) with this property. Many iterators have 



this property, such as Gauss-Seidel or weighted Jacobi (with weight w = 2/3 in ID or = 4/5 in 2D), but 



not Jacobi (see [41', pag. 30-32] for more details). From now on, by (17) we shall intend the Gauss-Seidel 
version of (0, (|8j), ([§, i.e.: 









1 




2 




= u, 



(m) 



+ /^AT At gb 



In order to explain the multigrid approach, we just describe the two-grid correction scheme (TGCS), because 
all the other schemes, such as F-cycle, IF-cycle, -F-cycle or Full multigrid cycle, can be easily derived from 
it (see |41t Sections 2.4, 2.6] for more details). The TGCS consists into the following algorithm: 

1. Set initial guess u/i = 

2. Relax v\ times on the finest grid: for k from 1 to v\ do 

3. Compute the defects 

Th = ih + A/iUfe 
gb = gb-gNi^'') 

4. Transfer the defect Th to a coarser grid with spatial step 2h by a suitable restriction operator 

^2h = Ik (r/^) 

5. Solve exactly the residual problem on the coarser grid 

- A2he2h = r2h (19) 
gfie2h) = ga (20) 

g^N{e2h) = h (21) 

in the unknown G S{Q.2h U T2h) 
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6. Transfer the error to the finest grid by a suitable interpolation operator 



7. Correct the fine-grid approximation 



8. Relax 1^2 times on the finest grid: for k from 1 to 1^2 do 



We have just to explain the steps concerning grid migration (steps 4 and 6). 
1.3 Transfer grid operators 

In this section, we describe the transfer grid operators for vertex-centered grid. We observe that our approach 
is based on the discretization of the equations on the various grids (both for inner and ghost points) . This 
approach is very different from algebraic multigrid. As a consequence, the interpolation and the restriction 
operators are not the transpose of each other. 

1.3.1 Restriction operator 

Since such operator will act on the defect G S{^h) (step 4), we must determine l2f^rh{x) for any x G Q2h 
using only values inside Qh- This is justified by the fact that the defect of the inside grid points (referred 
to the Poisson equation) may be very different (after few relaxations) from the defects ga, gb (referred to 
the boundary conditions and stored computationally in the ghost points), because the operators (for inner 
equations and for boundary conditions) scale with different powers of h. Then, let x G i^2h and refer to Fig. 
[2] (upper part). If x is not near an outside grid point, i.e. min{|x — a\ ,\x — b\} > h, then we will use the 
standard full- weighting restriction operator (FW): 




(22) 



while a x — h < a or X + h > b we set respectively 



Ik^hix) 



- {rhix) + rh{x + h)) 



(23) 



or 



- {rhix - h)+rhix)) . 



(24) 
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2h 



Restriction operator 

















1 — \ — \ — f- 


H 










— 1 


i — 1 — 1 


1 ^ — 

-1 

c 


T T T T 


V 1 



^1. 



^1 



2h ^ 



Interpolation operator 

^ — + — + — 1 — 'i — \ — t — t- 




H 1 



Fig. 2: Vertex-centered discretization in ID. Inner grid nodes (red circles) and ghost points (empty circles) on the fine and 
coarse mesh. The dashed lines represent the action of the restriction (up) and the interpolation (down) operators. 



1.3.2 Interpolation operator 

Since the interpolation operator acts on the error (step 6), which is continuous across the boundary, we do 
not need to separate the interpolation for inner equations from the interpolation of ghost points, and then 
we just use the standard linear interpolation operator (see the lower part of Fig. [2]): 

Iff'e2h{xj) = e2h{xj) if j is even 

Il''e2h{xj) = l{e2hixj-i) +e2h{xj+i)) if j is odd. 

Remark. 1 (y-cycle) The y-cycle algorithm is easily obtained from the TGCS recursively, namely 
applying the same algorithm to solve the residual equation in step 5. To terminate the recursion, an exact 
solver is used to solve the residual problem when the grid becomes too coarse. 

Remark. 2 (Ty-cycle) The VF-cycle is similar to the T^-cycle, with the only difference that the residual 
problem is solved recursively two times instead of one (in general schemes, 7 times, but 7 > 2 is considered 
useless for practical purpose). 

Remark. 3 (Coarser operator) We observe that the discrete operator A2h in step 5 is just the 
operator obtained discretizing directly the continuous operator in the coarser grid, and not the operator 
obtained by the Galerkin condition 

= hh h ■ 

The latter approach, typical of algebraic multigrid, makes the algebraic problem more expensive from a 
computational point of view and does not take advantage of the fact that the discrete problem comes from 
a continuous problem. 
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2 High-dimensional case 



In this case the defect of the boundary conditions has to be transferred in a suitable way to a coarser grid. 
The restriction has to be performed separately from the restriction of the interior equations, since these 
defects may show a sharp gradient crossing the boundary, because the discrete operators scale with different 
powers of h. 

In case of arbitrary domain, ghost points may have a complex structure and the restriction cannot be 
defined straightforwardly as in the rectangular case, where ghost points are aligned with the grid and the 
restriction can be performed by a one dimensional operator. 

For arbitrary domain we first need to extend the defect in a narrow band outside the domain constant 
along normal direction, and then we can operate the restriction as in the interior of the domain. For the sake 
of clarity, we describe the multigrid strategy in the two-dimensional case, but the procedure can be extended 
straightforwardly in more dimensions. We always refer to the second order method proposed in [13j, which 
is briefly recalled here. 

2.1 Model problem and relaxation scheme 

Let D = [—1, 1]^ be the computational domain, C -D be a domain such that n dD = 0. Let T£),Tn 

o o 

be a partition of 50 (i.e. F^) U F^v = dQ, F_d H Tn= where the interior points are computed in the d — 1 
dimensional topological space). Consider the model problem 

—An = / in Q 

u = qd on F^, (25) 
gN on Fat, 



du 
dn 



92 92 

where n is the outward unit normal, A = -—- + -—- is the Laplacian operator, / : — )■ M, (7/5 : F/j 
qn- Fat — )• M are assigned functions. 



In order to solve the elliptic problem (25), we can transform it in an evolutive problem (with a fictitious 

2t problem: 

An + / inn (26) 



time) that we call the associate time- dependent problem: 

du 



^ = fJ-oign - u) on Yd (27) 

du f du\ ^ , , 

f^N \ gN - TT ] on Fat (28) 



dt \ dn ^ 

u = uq in il, when t = 0, (29) 

where and /iTv are two positive constants. Then we look for the steady state solution. An iterative 
scheme can therefore be obtained by discretizing the associate time- dependent problem and considering the 
time just as an iterative parameter. 
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2.2 Level-set function 



In order to keep track of the boundary F, we introduce the level set function (/?o : -D — )• M, in such a way: 

o 

(x, y) en^^ (po{x, y) < 0, (x, y) G dn cpo{x, y) = 0. 



The outward unit normal to the boundary is 



n 



General references on the level set method for tracking interfaces are, for examples, [31] or [25j. From the 
level set function, we can obtain the signed distance function ip by fast marching methods [25] or by the 
reinitialization procedure based on the numerical solution of the following PDE 

^ = sgn((^o)(l-|V^|), (30) 

as we can see, for instance, in |40| [36l I17j . A signed distance function is preferred to a simple level-set 
function because sharp gradients are avoided and it is simpler to compute the boundary closest point to a 
given ghost point. Now we assume that |V</?| = 1 and suppose we know the signed distance function just at 



the grid nodes. In practice, Eq. (30) has to be solved for a few time steps, in order to compute the distance 



function a few grid points away from the boundary, z 



2.3 Relaxation operator 



Let us introduce some notation. Let d E N be the dimension of the problem, > 1 be an integer and 
h = 2/N the spatial step. Let Dh = = {ji, . . . ,jd) G {—N,N}'^ and Jl/i = O n Dh be the discrete 
versions of D and VL respectively. is the set of grid points. Two points x', x" in are called neighbor if 
I^j=i x'j ~ x'j = h. We call ghost point any grid point that is both outside Q. and neighbor to a grid point 
inside VL. We call Th the set of all ghost points. Let Ih be a general subset of Dh. We introduce the linear 
space of grid functions over and we denote it S[Ih) = {^h '■ Ih 

From now on, we shall consider d = 2, but the results are valid also for d > 2. 

Then, we write the basic iterative scheme {relaxation scheme) discretizing the time- dependent problem 



(26)-(29). For any grid point {jh,ih) of Qh, we write an equation obtained from the discretization of (26) in 
such point, using forward Euler in time and central difference in space and taking the maximum time step 
consented by the CFL condition, i.e. At = /i^/4 (in general it is At = h'^/{2d)): 



(m+l) _ -. /. f ,2 , ^(m) ,^,(m) , . ^(m) 



- 1/4 (h^f,, + + + <]u + <-u ) • (31) 



Eq. (31) is equivalent to discretize directly the first equation of (25) using central difference in space and 



applying Jacobi iteration scheme. 
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Since we have used the standard 5-point stencil even for grid point close to the boundary, we have added 
new unknowns to the linear system (ghost points). 

To close the system of equations (31), we must write one equation for each ghost point. This can be 
done in three simple steps. Let G = (xcUg) be a ghost point. 

1. Making use of the signed distance function cp, we can compute the closest boundary point to G, that 
we call B (see Figure Isl), by: 



B={xB,yB) = G-nG-ip{G) = G 



W7\ 



(32) 



G 



using a second order accurate discretization for V(/7, such as central difference in G. 
2. Compute the nine-point stencil (depicted in Fig. [3]) in Upwind direction, i.e.: 

Stc = ^{xG + Sx ki h,yG + Sy k2 h): {ki, k2) G {0, 1,2}^| , 

where = sgn(xB - xg) and Sy = sgn{yB - yc)- 



(33) 



3. Let CstcVA be the biquadratic interpolant of the numerical solution u in the stencil Stc- If -B G F^i, 



the iteration for the ghost point G will be obtained from the discretization of (27): 



u. 



(m+l) 
G 



(34) 



while if i? G Fjv, the iteration for the ghost point G will be obtained from the discretization of (28) 



u, 



(m+l) 
G 



(m) 
G 



B 



(35) 



The constants /^d and are chosen in order to satisfy a CFL condition, i.e. fio^t < 1 and /hat At < 
2/i/(3\/2) [13]. In numerical tests of Section^ we choose fiD^t = 0.9 and /XArAt = 0.9 • 2/i/(3\/2). bmce 
At = h'^/A, then = S.6/h'^ and fiN = 7.2/(3\/2/i). 

Remark 1 (Accuracy of (32)). The accuracy of the evolution of point B in (32) depends on the 
accuracy at which if is computed. If ip is known to order h^, p G {2,3}, then B will be computed to the 
same order of accuracy, provided we are far from singularities in ip. 

In Eq. (32) one could omit the term |V(/?| in the denominator, because |V(^| = 1 if is a signed distance 
function. However, it is better to keep such term, in case ip is only approximately a signed distance function. 

Remark 2 (Upwind stencil). The reason for which we use an Upwind stencil StG is simple. Let us 
rewrite the Neumann boundary condition (28) as: 



du 
dt 



+ ^Arn • Vu = IJ.N9N, 



(36) 
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Fig. 3: B is the boundary closest point to G, while the red points are the nine-point stencil in Upwind direction referred to the 
ghost point G. 



in such a way it appears to be an hyperbolic equation which propagates the solution u along the characteristic 
(normal direction to the boundary) with speed fiN t^- Then, it is preferred to use an Upwind stencil in 
direction — n to discretize the spatial term, in order to guarantee convergence (see [39] for more detail about 
Upwind schemes in Conservation laws). 

Remark 3 (Reduced stencil). If Stc is not fully contained in U Th, therefore we reduce the 
nine-point stencil to a smaller stencil, such as a 2 x 2 stencil or a (less accurate, more robust) three-point 
stencil. Such a reduction occurs rarely, and does not degrade the whole accuracy of the method (see [13j). 

Using the simplified notation, the iterative scheme converges to the solution of the problem: 



• Ufi £ S{Qh U Th) is the unknown; 

• Ah : S{Clh U Tfi) — )• S{Clh) is the standard discrete version of the Laplacian operator, namely: 

Ah^hix, y) = (w/i(a; + h,y)+ w/,(x - h,y) - 4w/,(x, y) + Wh{x, y + h) + Wh{x, y - h)) 

for any S ^(il/i U Th) and {x,y) G VLh] 

• ih€ S{^h) is defined by /h(-P) = /(P) for any grid point P G Qh] 




(37) 



where: 
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Lh'- S{Qh U Th) —7- S{Th) is the discrete version of boundary conditions, namely: 

if BgTd 



if Set 



TV 



(38) 



B 



for any w/j E S{^h U T/,) and G E T/^; 
• gft E 5(r;i) is defined by: 

for any ghost point G ^T^. 
2.4 Multigrid approach in 2D 



gD{B) if BeTo 
gN{B) if S E r^v 



Consider the Poisson problem ( 25 ) and suppose we have a discrete approximation of the form ( 37 ) . Therefore 



we are dealing with non-eliminated boundary conditions. Let us introduce, for any spatial step h, an exact 
solver 

vih = Sh{h,^h) (39) 



of the system (37), and denote by 

^h- si^h^Th) X s{^h) X ^(r;,) 

the relaxation operator, namely the iterative scheme 



u 



(m+l) 



(40) 



(41) 



(31), (34), (35). 



converges to the solution of (37) as n — )• +oo. In details, the iteration (41) summarize the iterative scheme 



As in one dimensional case, we will intend by (41) the Gauss-Seidel version of (31), (34), (35), in order 



to deal with a proper smoother, and we have to order all points of Q.h U in some way. Let us choose the 
lexicographic ordering (GS-LEX): 



{x',y')<{x",y") 



x' < x" or x' 



x", y' < y". 



The relaxation scheme can be easily extended to more efficient kinds of smoothers, such as Red-Black 
Gauss-Seidel (see [41j): however, we limit ourselves to study of the GS-LEX smoother. 

In order to explain the multigrid approach, we just describe the two-grid correction scheme (TGCS), 
because all the others schemes, such as ^-cycle, VK-cycle or Full multigrid, can be easily derived from it 
(see |4H Sections 2.4, 2.6] for more details). The TGCS consists into the following algorithm: 
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1. Set initial guess Uh = 



2. Relax vi times on the finest grid: for k from 1 to vi do 

3. Compute the defect r/j E S{i^h U T^) such as: 



rh 



fh + Ah Uh in ^Ih 
9h - Lh u/i on Th 



4. Transfer the defect to a coarser grid with spatial step 2h by a suitable restriction operator 

5. Solve exactly the residual problem in the coarser grid 

e2/i = S2h {^2h>'^2h) 

where r^^, = r2h\n^^ and rg = r2h\r^^; 

6. Transfer the error to the finest grid by a suitable interpolation operator 

7. Correct the fine-grid approximation 

Uh-. = Uh + Bh 

8. Relax 1^2 times on the finest grid: for k from 1 to 1^2 do 

Uh:=^h (u/i,f/i,g/,) 

We have just to explain the steps concerning grid migration (steps 4 and 6). All the other steps are clear. 
2.5 Transfer grid operators 



If (40) has the smoothing property, after ui relaxations (step 2 of the algorithm) we have a smooth defect 
r/i. Therefore, we can hope to transfer this defect to a coarser grid without losing much information. The 
defect Th defined in the step 3 belongs to S{Qh^^h)- In order to transfer it to a coarser grid, it is convenient 
to extend in some way this defect in the whole computational domain Dh (i.e. r/, G S{Dh)), in such a way 
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we can use the standard full- weighting stencil for the restriction operator /j^^: S{Dh) 
(see mi pag. 42]) 

"12 1 



^2h 



1 

16 



2 4 2 
1 2 1 



5(Z?2/i), that is 
(42) 



2h 



In general, by the stencil notation 



h,-i hfl til 



2h 



we will intend the restriction operator /2/i defined by: 

l2h^h{x,y)= ^ tijWhix + jh,y + ih), 

where only a finite number of coefficients tij is different from zero, and Rk = {—k, . . . , k}'^ for some positive 
integer k. In practice k = 1 allows second order restriction operator. 

Let us suppose we have extended the defect to the whole computational domain (as it is carefully 
described in Sec. 2.6). Anyhow, since we have different operators for inner equations and for boundary 
conditions, the defect is smooth separately inside Qh and along the ghost point (or — because of 
the extension) , but it is not smooth in all i^h U Th (it shows a sharp gradient crossing the boundary, as we 
can see in Fig. |4]). 

For this reason, it is convenient to transfer separately on the coarse grid the defect in and in Dh — ^h- 
To do that, we introduce a partial grid transfer 

4,: S{Dh) X V{Dh) X S{D2h) S{D2h), 
where V{Dh) is the family of all subset of D^. Roughly speaking, 

^2h = l2h {^hJh,W2h) 

means that we transfer w/j to a coarser grid using only the points of I/i, leaving unaltered the value in 
the points out of Ih and already stored in 'W2h (to better understand, we can think Ih = VLh). 

In details, let (x, y) G D2h H Ih- We focus our attention to the neighborhood of (x, y), that is M{x^ y) = 
{{x + jh,y + ih): j,i = -1,0,1}. 

Now consider the maximum full rectangle T with vertices belonging to AA(x, y) and such that TCi C 
M{x,y) n Ih (see Fig. [s]). Therefore, the stencil we use in {x,y) to transfer w/j to a coarse grid depends on 
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60 



70' ° 



Fig. 4: After few relaxations, the defect show a sharp gradient crossing the boundary (the defect is not yet extended outside), 
especially for the Neumann boundary. Here the domain is a circle. 



the size of T- In fact, let TO be a 3 x 3 points (i.e. J\f{x, y) C /^). In this case we can use the standard 
full- weighting stencil ( 42 ) . 

Now let T n Dh be a 3 X 2 points. Without loss of generality, we can suppose the vertices of T are 
(x + jh, y + ih), with j E {—1, 0}, i G {—1, !}• In this case, the stencil we will use is: 



{ik'^hj ix,y) = ^ 



2 2 
4 4 
2 2 



ix,y), 



(43) 



2h 



while, if T is a 2 X 2 points, with vertex {x + jh, y + ih), j, i EE {—1, 0}, the stencil will be: 

h 



i2h^h ] (x, y) 



1 

16 




4 4 
4 4 



{x,y), 



(44) 



2h 



This three case are summarized in Fig. [s] (where J/j = ^h)- 

Finally, for all points {x,y) E - we set ^2h{x,y) = W2h{x,y). 
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Fig. 5: Upper, the nine points of Mix, y) and the green boundary of the rectangle T. The bold red point is on the coarser and 
finer grids, while the little red points are on the finer grid. The arrows represent the action of the restriction operators. Below, 
the respective stencil to be used. 



In such a way, we can easily define the restriction operator '■ S{Dh) S{D2h) as foUows: 

l2h^h = Ik (^h, Dh - Qh, (wh, ^h, 0)) . 



(22), (23), (24). 



Note that the stencils (42), (44), (43) can be derived as tensor products of the one-dimensional restriction 



2.6 Extension of the defect 

In every ghost point we store the defect of the boundary condition concerning that ghost point. In formulas, 
we have seen in step 3 of the TGCS algorithm that rh(G) = (g/j — L^Uh) (G), for any ghost point G. But 
g/i(G) = g{B) and LhUh{G) is the reconstructed boundary condition in B of the boundary operator L (see 
(|38[)), where B is the closest boundary point to G (i.e. the orthogonal projection on the boundary). In 
summary, the defect is stored in a ghost point G, but it is geometrically referred to a boundary point B 
placed along the normal direction. When we switch to a coarse grid, some ghost point Gi may not be ghost 
point in the fine grid, i.e. T2h CI is not true (see Fig. [g]). 

Then, no acceptable value of the defect is stored in Gi. Indeed, we expect that r2h has in the ghost 
point Gi the defect of the boundary conditions referred to Bi. Hence, if we extend the defect rh outside 
Qh constant along the normal lines to the boundary, we will find rh (Gi) as an approximation of the defect 
of the boundary conditions in Bi. After coarsening (performed using only points outside as described 
before), the ghost points of the coarser grid will contain the expected values of the defect. 
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Fig. 6: Red bold and small points are grid points of Qh, while red bold points are grid points of Q2h- G\ is a ghost point on 
the coarser grid, but not on the finer grid, then no value of the defect is stored in it. and Qy are the two upwind near points 
to Gi. 



The extension of the defect is performed by solving the transport equation 

dr 



dr 



+ Vr • n = 



in a few steps of a fictitious time r, where n = (n^^, Uy) = V^j/I V(^| is the unit normal vector to the level-set, 
while r{x^ y) is a continuous version of r/j (i.e. r(x, y) is a continuous function defined va D — Vt and such that 
r[G) = Vh{G) for any ghost point G). In details, we compute few steps of the following iteration scheme: 



(m+l) 



(P)=rf)(P) + 



Ar 



(45) 



for all P G -D/j — (il^ U F/i), where Qx and Qy are the two upwind near points to P, i.e. 

Qx = P - sgn(n^) hi, Qy = P - sgn{ny) /i j. 



However, it is sufficient to perform the iteration (45 ) only in a narrow band with width 3 h. In order to speed 
up the extension, we can perform (45) in a Gauss-Seidel fashion, sorting points in — {^h U Th) by the 
distance from the boundary (it can be done using the distance function ip), and starting the computation 
r^f^~^^\p) in (45) from the closest ghost point P to F. 
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2.6.1 Interpolation 

Since the interpolation operator acts on the error, which is continuous across the boundary, we just use the 
standard bihnear interpolation operator: 



T2h _ 1 



1 2 1 

2 4 2 
1 2 1 



■2h 



3 Numerical tests 

In all the following numerical tests we always choose the Dirichlet and Neumann parts of T = dVt as: 

Yd = {(x,y) G r: X < 0}, Tn = {{x,y) G T: x > 0} . 

The Local Fourier Analysis (LFA) is a powerful tool to obtain the theoretically convergence factor by 
analyzing separately the action of different parts of the multigrid algorithm to high and low frequency 
components of the error. For a detailed explanation of the LFA, we refer to [41| Ch. 4]. 



Before to apply the LFA, one has to be sure that the relaxation operator (40) has the smoothing property. 
Roughly speaking, the smoothing property is the property to dump high frequency components of the error, 
in order to make it smooth after few relaxation sweeps. 

When the multigrid algorithm applies to a regular rectangular domain, the LFA and smoothing analysis 
are well studied. In the case of arbitrary domain, as Achi Brandt points out in |41l pag. 587], there are 
some boundary related difficulties about the discretization and relaxation near the boundary: 

• There is no a general smoothing analysis when the boundary is not aligned with the grid; 

• The residuals should be reduced near boundaries more than in the interior; 

• The coarsest grid has not to be too coarse, because it should catch the curvature of the boundary in 
order to guarantee the convergence. 

Now, we perform a numerical test in order to check if the convergence factor is close to the predicted one by 
LFA, which is obtained for rectangular domain with periodic boundary conditions, i.e. without taking into 
account boundary effects. Note that the multigrid algorithm described before may be seen as an iterative 
scheme: 

(m+l) _ , . (m) , 

for some matrix and vector b/j. We call p the convergence factor, which is the spectral radius of 
the matrix M^. For rectangular domain with periodic boundary conditions and constant coefficients, the 
convergence factor is said to be local and it is denoted by pioc- The convergence factors predicted by LFA 
for Gauss-Seidel LEX relaxation and FW restriction operator are listed in Table [l] (see [HI pag. 117]). 
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Table 1: Predicted convergence factor pioc by LFA for GS-LEX and FW restriction operator. 



U = Ui + 1/2 


1 


2 


3 


4 


Ploc 


0.400 


0.193 


0.119 


0.084 



In all the numerical tests we perform, the convergence factor is estimated as the ratio of consecutive 
defects, i.e.: 







oo 




oo 



p = p 



for m very large. In order to avoid difficulties related to numerical instability related to the machine precision, 
we will always use the homogeneous model problem as a test, namely (25) with f = Qd = dN = 0, and 
perform the multigrid algorithm starting from an initial guess different from zero. Since we are just interested 
at the convergence factor and not at the numerical solution itself (which approaches zero indefinitely for 
homogeneous problem), a reasonable stopping criterion will be 

Ipim) _ „(m-l)| 

' ^ < 10^1 



Note that, since we want to study the efficiency of the multigrid components proposed in this paper 
(smoother, restriction, ...), it is sufficient to study basic kind of multigrid such as V-cycle and W-cycle, 
while a more efficient algorithm (such as FMG) can be easily derived. 

3.1 ID numerical test 

Referring to Sec. [l| let us choose [a,b] = [—0.743,0.843] C [—1, 1]. The finest grid is obtained dividing the 
whole computational domain [—1, 1] into = 64 subintervals, while the coarsest grid is obtained dividing 
[— 1, 1] into Nc = 8 subintervals. The computed convergence factors are very close to those ones predicted 
by LFA (Table [T]), namely p = 0.185 for i/ = 2 and p = 0.122 for u = 3. 



3.2 An initial test in 2D 



We start testing the multigrid algorithm on a circular domain Q with center (\/2/20, \/3/30) and radius 
r = 0.563 (Fig. [7]). 

The measured convergence factors for TGCS, V^-cycle and IF-cycle are listed in Table [2| 
As we can see, the measured convergence factor is far from the predicted one by LFA (Table [l]). Then, 
some boundary effect degrades the convergence factor. Note that in ID such boundary effects do not degrade 



the convergence factor (Ex. 3.1), because we have only two boundary points, and the degradation is due to 
the oscillating behavior of the residual on the tangential direction to the boundary, that does not exist in 
ID. Then, in 2D we must smooth the error also along the tangential direction to the boundary. To overcome 
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Fig. 7: Circular domain and the coarsest grid used to capture the curvature. Blue ghost points refer 
to Dirichlet condition, while green ghost points refer to Neumann condition. Red lines are normal to 
the boundary. 



Table 2: Measured convergence factor p with z/i = z/2 = 1 on the left and with i/i = 2, 1^2 = 1 on 
the right. 



N 


TGCS 


F- cycle 


VF-cycle 


64 


0.67 


0.68 


0.71 


128 


0.68 


0.73 


0.68 


256 


0.70 


0.71 


0.70 



N 


TGCS 


F-cycle 


W-cjcle 


64 


0.58 


0.72 


0.58 


128 


0.58 


0.73 


0.59 


256 


0.61 


0.83 


0.60 



this difficulty, we apply, after a single relaxation and at each grid level, A extra relaxation sweeps on all 
ghost points Th and on all inside grid points of ilh within 5 > distance from the boundary (the extra 
computational work is 0{N), then negligible as A'^ — t- 00). It can be proved numerically that a good choice 
of these parameters will be: 

A = 5, 5 = 3h. 

The explanation of the optimal value A = 5 is the following: the degradation observed in Table [2] is an 
indication that the error decays much slower at the boundary. Assuming that the convergence factor in 
Table [2] is essentially the convergence factor at the boundary, p^, we want to match it with the convergence 
factor at the bulk, therefore Xopt is the smallest value of A for which < pi- The value pi, in turn, can 
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be computed as the convergence factor for large value of A. 

Investigating the smoothing property, we observe that choosing the initial error as an high frequency 
component, the error is not smoothed after few relaxation sweeps. While, if we add the extra-relaxations, 
the error become sufficiently smooth (Figs. [8]j9|. 

3.3 Some numerical results 

In this section we confirm numerically the improvement of the convergence factor if we apply extra- 
relaxations, and we compare the relaxations with other well-knowns alternative such as the Kaczmarz and 
the block relaxations. In all numerical tests, we choose an arbitrary domain O assigning a level-set function 
(fQ. Then we reinitialize it by the procedure described in Section [2. 2[ obtaining the signed distance function 
Lp. Afterwards, we perform the multigrid technique applying the M^-cycle algorithm instead of the V-cycle, 
to ensure the independence of the convergence factor p from the step size h (as explained for example in |41|, 
pag. 78]). Several tests are performed for each domain, based on the different size of the finest and coarsest 
grids. The finest grid is obtained dividing the whole computational domain D into subintervals in each 
Cartesian direction, while the coarsest grid is obtained replacing N with Nc- The solution on the coarsest 
grid is obtained by a direct solver. 



3.3.1 Circular domain 



In this case we can choose as a level-set function directly the signed distance function, which is known 
analytically: 

ip{x, y) = \J{x- \/2/20)2 + {y- V^/mf - 0.563. 

The zero level-set is represented in Fig 10 (top- left). Different value of the convergence factor are listed in 
Table [3] (for i/ = z/i + f 2 = 2 and = 3). They are really improved with respect to those obtained without 
extra-relaxations (Table [2]) for the same test. 



3.3.1 



We use N X N number of grid 



Table 3: Convergence factor of the numerical test of Section 
points in the finest grid; Nc x Nc number of grid points in the coarsest grid. Left table: v — 1^1+1^2 = 2, 
right table 1/ = ui + 1^2 — 3. 



N 

Nc 


16 


32 


64 


128 


256 




N 

Nc 


16 


32 


64 


128 


256 


8 


0.052 


0.053 


0.11 


0.13 


0.14 




8 


0.06 


0.03 


0.09 


0.08 


0.08 


16 




0.061 


0.11 


0.13 


0.14 




16 




0.04 


0.09 


0.08 


0.08 


32 






0.11 


0.13 


0.14 




32 






0.09 


0.08 


0.08 


64 








0.13 


0.14 




64 








0.09 


0.08 


128 










0.14 




128 










0.09 
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Fig. 8: High frequency initial error after 1 (up-left), 3 (up-right), 5 (down-left), 10 (down-right) relaxation sweeps and A = 
extra-relaxations. 




Fig. 9: High frequency initial error after 1 (up-left), 3 (up-right), 5 (down-left), 10 (down-right) relaxation sweeps and A = 5 
extra-relaxations. 
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3.3.2 Comparison with the Kaczmarz and the block relaxations 

Note that the relaxation scheme ([7|), Q, ^ is composed by a Gauss-Seidel iteration over inner grid points 
and a suitable relaxation over ghost points (boundary conditions). As an alternative to the relaxation 
of the boundary condition, we can use the Kaczmarz relaxation [26] near the boundary, which is known 
to be unconditionally convergent. Let us recall the Kaczmarz iteration scheme for a subset of equations 
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J7 C { . . . , A'i + Ng} of a linear system Lu = f: 



for j £ do: u 



TEMP . 



TEMP 



U 



+ 



/J 



TEMP 



'"I ) 



'■J\\2 



The symbol < •, • > denotes the inner product operator and Ij is the j-ih row of the matrix L. If we 
choose J = {. . . ,Ni + Ng} then we obtain the classical Kaczmarz relaxation scheme for the solution of the 
linear system Lu = f, and the iteration scheme is equivalent to a Gauss-Seidel relaxation for the system 
L^Lu = L^f. In our case, one iteration of the alternative relaxation we want to study is composed as 
follows: we perform a Gauss-Seidel sweep in the interior of the domain, followed by A Kaczmarz iterations 
over ghost points and inner points close to the boundary (say within S distance from the boundary). 

Another alternative is represented by the block relaxation [16j. As we point out in [T3], the elimination 
of the boundary conditions is hard to perform in high dimensions, while in one dimension it is a trivial task 
and leads to a diagonally dominant linear system. A middle ground between the elimination of the boundary 
conditions and the relaxation operator we use in this paper is the block relaxation. Let us describe it in 
details. For each grid point P G Jl/j U F/i we choose a stencil Stp C U T^. For instance, ii P £ we 
choose Stp = StpQ (1 (Jl/i U F/i), where Stpg is the 3x3 stencil centered at P, else if P € F/i we choose the 



stencil Stp defined in (33). One iteration of the alternative relaxation is composed as follows. We perform a 



Gauss-Seidel sweep in the interior of the domain except in grid points within 6 distance from the boundary. 
For each grid point P £ QhUTh within 6 distance from the boundary we rewrite the linear system Lu = f 
as follows (by a permutation of rows): 



Lu = f <^ 



^1,1 

A2A 



^1,2 
^2,2 



Ul 
U2 



fl 
f2 



where ui is referred to those grid points belonging to Stp. Therefore, we update the values of ui as: 

Ul = A^^ {fl - Ai^2U2) ■ 

We perform a comparison between the relaxation proposed in this paper (that we call new iteration in 
the following plots) and the two alternative relaxation described above. Such a comparison is carried out 
in terms of smoothing factor and convergence factor. We perform the comparison using the TGCS for the 
test case of the circular domain 13.3.11 with N = 64. 



In Fig. 11 we plot the smoothing factor /x for the three iteration schemes, which is estimated by the ratio 
of subsequent defects after each iteration, i.e. 



(m) 



„(™-i) 
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In practice, we perform only the iteration schemes, without taking into account the effects of the multigrid 
procedure. In order to better capture the behavior of the smoothing factor, we choose an initial guess being 
highly oscillant, for example u = sin(407rx) sin(507ry). 

In Fig. 12 we depict the convergence factor p for the Kaczmarz and the new iteration against the 



number of extra-relaxations A (for comparison, we also plot the convergence factor of the block relaxation as 
an horizontal line, since it does not depend on A). After five extra-relaxations, the new iteration reaches a 
plate configuration, since it achieves the convergence factor of the Gauss-Seidel smoother for inner equations, 
i.e. the convergence factor predicted by the LFA (see Table [T]). The Kaczmarz iteration falls down slower, 
while the block iteration already provides the optimal convergence factor. The computational cost of five 
point-iterations of the new method is considerably lower than the cost of one block-iteration. 



1 

0.9 

0.8 - 

0.7 - 

0.6 - 
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- 



□ Kaczmarz 
* Box relaxation 
° new iteration 



ja^-o-o- 




5 10 15 

m-th iteration 

Fig. 11: Smoothing factor /i against the 
number of iterations for the three iterative 
schemes: Kaczmarz relaxation, Block relax- 
ation, new iteration. 



20 



1 

0.9 
0.8 
0.7 
0.6 
0.5 
0.4 
0.3 
0.2 
0.1 



□ Kaczmarz 

Box relaxation 

° new iteration 



t3- 



^-e — e — e — 9 — & — s - - o — 9 — e — g- 



10 15 



Fig. 12: Convergence factor p of the entire multi- 
grid against the number of extra-relaxations A for 
the three iterative schemes: Kaczmarz relaxation. 
Block relaxation, new iteration. 



3.3.3 Ellipsoidal domain 

The level-set function is: 



_ {X{x, y) - ^/2/20)2 (y(x, y) - ^3/30)^ 
0.5632 + 0.2632 
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where 



y) = cos(7r/6) x — sin(7r/6) y, X{x, y) = sin(7r/6) x + cos(7r/6) y, 



and the zero level-set is represented in Fig. [To] (top-right). The convergence factor obtained are listed in 
Table [4] (for v = vi + V2 = and v = ?>). We observe as the convergence factor degrade choosing a coarsest 
grid too much coarse, but starting from a certain level of coarsest grid it is relatively close to the predicted 
convergence factor by LFA (Table [l]). 



Table 4: Conver gence factor of the numerical test of Section |3.3.3| We use N x N number of grid 
points in tlie finest grid; x number of grid points in the coarsest grid. Left table: = i/i +1/2 = 2, 
right table 1/ = ui + — 3. 



N 

Nc 


16 


32 


64 


128 


256 




N 

Nc 


16 


32 


64 


128 


256 


8 


0.34 


0.09 


0.14 


0.14 


0.15 




8 


0.44 


0.06 


0.12 


0.11 


0.09 


16 




0.65 


0.45 


0.19 


0.15 




16 




0.55 


0.30 


0.09 


0.09 


32 






0.14 


0.14 


0.15 




32 






0.13 


0.10 


0.09 


64 








0.15 


0.15 




64 








0.12 


0.08 


128 










0.15 




128 










0.09 



3.3.4 Saddle-shaped domain 

The level-set function is: 



ip{x,y) = 9 



1 



V3 



+ 



3\/3 3 



sm 



SVS 3 



and the zero level-set is represented in Fig. 10 (bottom- left). The convergence factor obtained for z/ 



z^i + 2^2 = 3 are listed in Table [s] (left). Also in this case, only if we choose = 16 and Nc = 8 (which 
actually is TGCS) the convergence factor is degraded. 



3.3.5 Flower-shaped domain 

The level-set function is: 

ip = r ■ 



2. ,3 



0.5 



+ 5x^y - lOx^y 



+ y2 



and the zero level-set is represented in Fig. 10 (bottom-right). The convergence factor obtained for v = 
v\^V2 = 'i are listed in Table [5] (right). This is the hardest numerical test, because of the indentation of 
the boundary. We need to start from a coarsest level = 32 to correctly capture the boundary profile and 
to make the discretization accurate. 
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Table 5: Convergence factor of the numerical test of Sections 3.3.4 (left) and 3.3.5 (right). We use 
N X N number of grid points in the finest grid; Nc x number of grid points in the coarsest grid. 
In this test we use u = vi + U2 = 3. For the flower-shaped domain, if A^c = 8 the multigrid does not 
converge. 



N 

Nc 


16 


32 


64 


128 


256 




N 

Nc 


16 


32 


64 


128 


256 


8 


0.36 


0.08 


0.09 


0.12 


0.09 




8 


n.c. 


n.c. 


n.c. 


n.c. 


n.c. 


16 




0.12 


0.09 


0.12 


0.09 




16 




0.89 


0.75 


0.50 


0.25 


32 






0.09 


0.12 


0.09 




32 






0.49 


0.25 


0.12 


64 








0.13 


0.09 




64 








0.24 


0.11 


128 










0.09 




128 










0.09 



Conclusion 

A multigrid technique for Poisson equation on an arbitrary domain and mixed boundary conditions is 
presented. This multigrid strategy can be applied to a general framework of ghost-point method in a 
regular Cartesian grid, in case of non-eliminated boundary conditions. Suitable transfer operators for inside 
equations and boundary conditions are provided. The convergence rate is improved by adding some extra- 
relaxations on the ghost points and in a narrow band of inside grid points close to the boundary. Numerical 
tests on different geometries have been performed. On simple domains (such as circle or ellipse) the optimal 
convergence factor is reached even if the method is used on very coarse grids, while for more complex domains 
(such as the flower-shaped one) the optimal convergence factor is obtained only on sufficiently fine grids. 
A comparison with other treatments of the boundary condition smoothing procedure has been carried out 
(Kaczmarz and Block relaxation), confirming that the smoother proposed in this paper is better in terms 
of convergence factor, and not worse in terms of smoothing factor. 

The application of this method to solve the pressure equation coming from the projection method of 
Chorin [TT] in the framework of the incompressible Navier-Stokes equation is in preparation. Several 
extensions of the discretization technique and multigrid approach are presently under investigation. We 
mention the case of discontinuous coefficients, which models, for example, a system composed by different 
materials separated by an interface; in such a case the method is suitably modified in order to achieve second 
order accuracy and a convergence factor being independent on the jump in the coefficient. A preliminary 
result can be found in [TJj. Another extension concerns the convection-diffusion equation in a moving 
domain, in order to study applications modeled by a Stefan-type problem. All these extensions will be 
coupled with the use of Adaptive Mesh Refinement to obtain accurate solution in the case of domain with 
complex boundary. 
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